
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%                                                                         %
%                           Direct Parameters                             %
%                                                                         %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%


% ---------------------------- Description -------------------------------- 

    % Given a vector of parameters, we can form a set of two equations,
    % each guaranteeing that one of the targeted moments is achieved,
    % depending on two variables: (i) alpha, (ii) I.

    % We use Matlab's fsolve to solve this non-linear system of equations.
    % As a result, two moments will equal their targets exactly: (i)
    % productivity growth rate, and (ii) aggregate entry rate.

%----------------------------- Preliminaries ------------------------------

    % Determine an initial guess for fsolve
    x0 = [0.6, 0.025];

    % Transforms the system of equations into a function of x only
    System = @(x) System_Functions(x,p);

    % Ask Matlab not to show any output
    options = optimset('Display','off');

    % Solve the system
    x = fsolve(System, x0, options);

    % Attributes the solution of the system to the appropriate variables
    p.alpha = x(1);
    p.In = x(2);

    % Clears the output of the system
    clear System options x x0;

%--------------------------------------------------------------------------
%                            Setting the System                           %
%--------------------------------------------------------------------------

    % Create an auxiliary function to build the system
    function F = System_Functions(x, p)

        % The first entry of the vector x: alpha
        a = x(1);

        % The second entry of the vector x: Innovation rate
        I = x(2);

        % ----------------------- Preliminaries ---------------------------

        % Based on alpha and delta, compute tau - see Proposition (M-3)
        tau = a*(p.eta+p.d0)/(1-a);

        % Based on alpha and wbar, compute qbar (to the sigma-1 power) - see Equation (M-5)
        qbar = a*p.l^(p.s-1) + (1-a)*p.wbar;

        % Determine the net rate of product accumulation - see Section (RP-2.2.2)
        psi = p.xstar - (tau + p.d0);

        % Determine the flow rate of entry z - see Proposition (M-3)
        z = (p.eta + p.d0)/(1 - a) - p.xstar;

        % ---------------------- First Equation ---------------------------

        % Compute the distance between productivity growth rate and target - see Section (RP-2.2.3)
        g_diff = p.eta/(p.s-1) + I + tau/(a*(p.s-1))*(qbar-1) - p.g_Moment;

        % ---------------------- Second Equation ---------------------------

        % Computes the numerator of Firm Age Density - see Equation (RP-4)
        Firms_function = @(y) psi*exp((psi-p.eta).*y)./(psi - p.xstar*(1-exp(psi.*y)));

        % Integrates the numerator over firm age for denominator - see Equation (RP-4)
        Firms = integral(Firms_function, 0, Inf);

        % Compute the distance between model and data Entry Rates  - see Section (RP-2.2.3)
        Ent_diff = z*((p.eta-psi)*Firms)^(-1) - p.EntryRate_Moment;

        % ------------------------- Stacking ------------------------------

        % Stacks up all equations to form a single vector - see Section (RP-2.2.3)
        F = [g_diff, Ent_diff];

end